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Polynomial model inversion control: 
numerical tests and applications 

C. Novara 


Abstract —A novel control design approach for general nonlinear 
systems is described in this paper. The approach is based on 
the identification of a polynomial model of the system to control 
and on the on-line inversion of this model. Extensive simulations 
are carried out to test the numerical efficiency of the approach. 
Numerical examples of applicative interest are presented, con¬ 
cerned with control of the Duffing oscillator, control of a robot 
manipulator and insulin regulation in a type 1 diabetic patient. 


I. Introduction 

Consider a nonlinear discrete-time system in regression form: 

yt = h{ut ,yt ( 1 ) 

ttt ~ (ttt—1, • ■ • , Ut—n) 

yi = {yt-i,---,yt-n) 

Ct — ■ ■ ■ j^t—n) 

where ut G U C M"" is the known input, yt G K"" is the 
measured output, G S C is an unmeasured disturbance; 
n is the system order; U and 2 = {^ S : ||^|| < f} are 
compact sets; the function h is Lipschitz continuous on 17^ = 
yn ^ jjn ^ where F is a compact set. U accounts for 
possible constraints on ut- 

Suppose that the system Q is unknown, but a set of noise- 
corrupted measurements is available: 

-D = {yuU^}t^.L ( 2 ) 

where the tilde is used to denote the samples of the data set 
V. 

Let 3^° C F" be a set of initial conditions of interest, TZ = 
{r = (ri,r 2 , ■■ ■)'. rt G F, Vf} a set of output sequences of 
interest, and 2" = {^ = (^i, ^ 2 , ■■■)'■ G S, Vf} the set of all 
possible disturbance sequences. 

The problem is to design a controller for the system 0 such 
that, for any ^ G 2, and for any initial 

condition yo G y^, the output sequence y = {yi,y 2 , ■ ■ ■) 
of the controlled system tracks any reference sequence r = 
(ri,r2,...) e TZ. 

To solve this problem, a novel data-driven control approach 
will be described in the following, based on the identification 
of a polynomial prediction model and on the online inversion 
of this model via the efficient solution of suitable optimization 
problems. A simplified version of the approach is presented 
in ||2l. 
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II. Data-based prediction model 
A model is considered, of the form 

y+ = /(u+,gr-) (3) 

where y^ = y^ is a prediction of the system output (over 
some finite time horizon), is a vector with the 

present and future input values and q~ = qf = {uf ,yf). The 
subscript indicating the time will be omitted in the reminder of 
the paper when not necessary. A parametric structure is taken 
for the vector-valued function /. In particular, each component 
fj of / is parametrized as 

N 

/f (•) = (•) (4) 

i=l 

where fi are polynomial basis functions, ofy are parameters to 
be identified and j = 1,..., rn^. The parameters can be 
identified from the data (|^ by means of convex optimization. 

III. Polynomial inversion control 

The proposed control approach is based on the on-line inver¬ 
sion of the model 0: at each time t > 0, given a reference 
sequence and the current regressor q~, a command se¬ 
quence M+ is looked for, such that the model output y^ is 
“close” to r+: 

y+ = f {u+,q~) =r~^. (5) 

Such a command sequence is found solving the optimization 
problem 

u* = arg min J (u, r~^,q~) (6) 

ueu^ ^ ' 

where 

J (u, r+, q-) = ||r+ - / (u, q~) II 2 + M HuH^ (7) 

and p, > 0 is a design parameter, determining the trade-off 
between tracking precision and command activity. 

The problem (|^ is solved at each sampling time, resulting in 
the following control law: 

< = ( 8 ) 

where is the first entry of the vector u* in (|^. 

The objective function Q is in general non-convex. Moreover, 
the optimization problem (|^ has to be solved on-line, and 
this may require a long time compared to the sampling time 
used in the application of interest. To overcome these relevant 
problems, three algorithms have been developed, allowing an 
efficient computation of the optimal command input for the 
following cases: 
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1) SIMO system and piecewise constant command input; 
the optimal solution can be computed “almost analyti¬ 
cally”. 

2) MIMO system affine in the cost function is convex, 
implying that the optimal solution can be obtained with 
“low” computational cost. 

3) General MIMO system, we will show below by means 
of extensive simulations that the algorithm is able to hnd 
always a solution “very close” to a global one, in very 
short times. 

The algorithms are based on a coordinate minimization ap¬ 
proach but are not described here. 


IV. Optimization algorithm performance 

EVALUATION 


The optimization problem]^ was considered, where /(•) is a 
polynomial function of degree dp and u G U C K"*, r^,q~ G 
K”*. This problem is analogous to (|^ but the dependence on 
time is not evidenced. The value /r = 0 was taken since, with 
this value, if is in the range of / (•), we know the global 
minimum of J (u, r+, g“) to be 0. 

Values of m in the set {1,2,4,6, 8} and values of dp in the set 
{1, 2,4, 6} were considered, corresponding to MIMO systems 
with up to 8 command inputs and models with polynomial 
degree up to 6. Note that in all the applications presented 
below, degrees 2 -F 4 led to a very satisfactory prediction and 
control performance. Degrees larger than 4 seem in general to 
not give any advantage. 

For each combination of m and dp in these sets, a Monte Carlo 
simulation was carried out, consisting of 50 main trials, each 
consisting of 100 sub trials (total number of trials: 5 * 4 * 50 * 
100 = 100 000 ). 

In each main trial, / (•) was dehned as a polynomial function 
of degree dp with sparse random coefficients. In particular, 
a number Ug of nonzero coefficients was assumed, with Ug 
ranging in the interval [0, 500] in function of m and dp (the 
nonzero coefficients were chosen according to a Gaussian 
distribution with zero mean and unitary variance). In each 
sub trial, a sequence was generated, 

where q~ and are vectors with random entries (chosen 

according to a uniform distribution with support [—1,1]), and 
i = 1,, 100. Then, for each i, the optimization problem (|^ 
was solved. Note that the decision variable u is different from 


the 


“true” input 


For each combination of the dimension m and the polynomial 
degree dp, the following indexes were considered to evaluate 
the algorithm performance: 

• -^2 = 5^ Ei° 1° {J (<> rf ,q~)-J (uf, r+, g")), 
where u* is the solution of the optimization problem 

computed for each random sample. Note that, in the 
present case, we know that J = 0. 

- J , rf ,q~)). 

2=i,...,5UUU 

• Tgc =average time taken by a Matlab .m function to solve 
a single optimization problem on a laptop with an i7 3Ghz 
processor and 16 MB RAM. The average was computed 
over the 5000 samples of the Monte Carlo simulation. 


• Tm ^average time taken by a compiled Simulink mex 
function to solve a single optimization problem on the 
same laptop. This function was generated in 10 of the 
50 main trials, since this operation is relatively complex. 
The average was thus computed over 10 * 100 = 1000 
samples of the Monte Carlo simulation. 

The obtained results are summarized in Table U It can be 
concluded that the coordinate descent minimization approach 
is able to hnd precise solutions (i.e., giving small values of 
the objective function) in short times for all the considered 
input dimensions and polynomial degrees. It can also be ob¬ 
served that using compiled mex functions allows a signihcant 
reduction of the computation times for problems involving 
polynomials with a not too high degree in u. A possible 
interpretation is that the Simulink automatic compiler looses 
efficiency for large degree polynomials. 



dp 


^2 

Eoo 

Tsc [s] 

Tm [s] 

1 

1 

3 

1.2e-14 

3.2e-14 

2.7e-4 

<1.0e-4 

2 

6 

1.9e-13 

1.9e-12 

3.0e-4 

<1.0e-4 

4 

15 

2.1e-13 

1.4e-12 

3.4e-4 

<1.0e-4 

6 

28 

l.le-13 

6.5e-13 

3.7e-4 

<1.0e-4 

2 

1 

5 

4.3e-12 

1.6e-ll 

8.7e-4 

<1.0e-4 

2 

15 

5.0e-3 

0.048 

1.6e-3 

1.4e-4 

4 

45 

4.1e-3 

0.022 

2.2e-3 

5.6e-4 

6 

81 

5.2e-3 

0.034 

5.4e-3 


4 

1 

9 

7.5e-5 

4.2e-4 

7.8e-4 

<1.0e-4 

2 

45 

8.2e-3 

0.039 

3.1e-3 

4.5e-4 

4 

116 

0.013 

0.047 

0.038 


6 

197 

0.014 

0.046 

0.17 

>Tsa 

6 

1 

13 

4.3e-4 

9.7e-4 

1.9e-3 

<1.0e-4 

2 

81 

0.013 

0.048 

0.011 

1.2e-3 

4 

197 

0.016 

0.048 

0.21 

>Tsc 

6 

339 

0.021 

0.049 

0.76 


8 

1 

17 

5.0e-4 

8.6e-4 

2.4e-3 

<1.0e-4 

2 

116 

0.019 

0.048 

0.10 

3.1e-3 

4 

289 

0.027 

0.049 

1.6 

>Tsc 

6 

500 

0.032 

0.049 

8.3 

>Tsc 


Table I 

Monte Carlo simulation results. 


V. Applications 

A. Duffing oscillator 

The Duffing system is a second-order damped oscillator with 
nonlinear spring, described by the following differential equa¬ 
tions: 

Xi = X2 

±2 = —aixi — a2x\ — 13x2 + u ( 9 ) 

y = xi+^ 

where x = {xi,X 2 ) is the system state (ti and X 2 are the 
oscillator position and velocity, respectively), u is the input, y 
is the output, and ^ is a zero-mean Gaussian noise having a 
noise-to-signal standard deviation ratio of 0.03. The following 
values of the parameters have been considered: ai = — 1, 
02 = 1, /3 = 0.2. For these parameter values and for 
certain choices of the input signal, this system exhibits a 
chaotic behavior, and this makes control design a particularly 
challenging problem. 

A simulation of the Duffing system (|^ having duration 400 
s was performed, using the input signal u{t) = 0.3sin(r) -I- 
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^“(r), wherer here denotes the continuous time and is 

a white Gaussian noise with zero mean and standard deviation 

0 . 2 . 

A set of L = 4000 data were collected from this simulation 
with a sampling period Tg = 0.1 s: 

yt}i=-1999 

where Ut = u{Tgt) are the measurements of the input and 
Vt = y{Tst) are the measurements of the output. 

A nonlinear controller was designed following the approach 
described in Sections [l^and [HI] This controller was applied 
to the Duffing system 

A testing simulation of the controlled system with duration 800 
s was performed, using zero initial conditions and a reference 
signal Tt generated as a sequence of random steps, filtered by a 
second-order filter with a cutoff frequency of 2 rad/s (this filter 
has been inserted in order to ensure not too abrupt variations). 
A Gaussian noise affecting the output measurements, having 
zero-mean and a noise-to-signal standard deviation ratio of 
0.03 was included in the simulation. In Figure [1] the output 
of the controlled system is compared to the reference. 

Then, a Monte Carlo simulation was carried out, where 
this data-generation-control-design-and-testing procedure was 
repeated 100 times. For each trial, the tracking performance 
was evaluated by means of the Root Mean Square tracking 
error 



The average RMS error obtained in the Monte Carlo simula¬ 
tion is HMS = 0.015. 

A simulation of the closed-loop system was also performed 
where vt = 0, Vf and was a step disturbance of amplitude 
0.5. The output signals obtained in these simulations are shown 
in Figure 

From these results, it can be concluded that the designed 
controller is able to (1) ensure a very accurate tracking, even 
in the presence of quite significant measurement noises; (2) 
reject/attenuate strong step disturbances. 




0 5 10 15 20 25 30 35 40 45 50 


Time [s] 


Figure 2. Above: disturbance rejection of the controlled system. Below: 
same figure, with zoomed y axis. Continuous (black) line: reference. Dashed 
(red) line: actual output. 


B. Robot manipulator 

The 2-DOF (2-degrees of freedom) robot manipulator depicted 
in Figure has been considered, where and (^2 the 
angular positions of the two segments of the robot arm, ui 
and U 2 are the control torques acting on these segments, li and 
I 2 are the segment lengths, and Mi and M 2 are the segment 
masses. The parameter values li = 0.8 m, I 2 = 0.7 m. Mi = 
2.5 Kg, M 2 = 2 Kg have been assumed. 

This robot manipulator is a MIMO system (with 2 inputs and 
2 outputs), described by the following continuous-time state- 
space nonlinear equations; 


i(r) 

y{T) 


= A‘=(z(t))z(t) -f B%z{t))u{t) 
zi{t) 

. ^2(t) _ 


( 10 ) 


where r is the continuous time, z(r) = [Ci('r) C 2 (''') Ci(t) 
(( 2 ]^ is the state, u(t) = W 2 (t)]~'^ is the input, and the 

expressions of A‘^(z(t)) G and B‘^{z{t)) G can 

be found in a. 



Figure 3. Robot Manipulator. 

A set of L = 5000 data was generated by simulation of ([T0|: 


Figure 1. Tracking performance of the controlled system. Continuous (black) 
line: reference. Dashed (red) line: actual output. 


'D = {yt,ut,}l= 


-4999 • 


The data were collected with a sampling time Tg = 0.02 s. 
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using the following input signals; 


= 


—20zj(T), if \zj{T)\ > 1.75 rad 
0, \fl<T<l + 500, I = 500,1500,2500,3500, 
and \zj{T)\ < 1.75 

U sin(wjiT) + U sin(wj2'r), otherwise, 

( 11 ) 

where j = 1,2, U = rand[50,150] Nm, uju = 

rand[0.05,0.09] rad/s, a;i2 = rand[0.5,0.11] rad/s a;2i = 
rand[0.04,0.1] rad/s wn = rand[0.7,1.2] rad/s. The notation 
U = rand[50,150] means that [/ is a number, randomly 
chosen according to a uniform distribution in the interval 
[50,150]. The feedback input on the hrst line of (111 was 
applied in order to limit the working range of zi and Z2 to 
the interval [—tt, tt] rad (the gain —20 and the threshold 1.75 
rad were chosen thorough several preliminary simulations). 
Measurement noises were added to yj, j = 1,2, simulated as 
uniform noises with amplitude 0.02 rad. 


From these data, two controllers were designed following the 
approach described in Sections|I^and|^ The hrst one is based 
on a general nonlinear prediction model. The second one is 
based on a prediction model affine in it+. For comparison, the 
controller in H has been considered, designed by means of a 
two-step method, consisting in LPV model identihcation and 
Gain Scheduling (GS) design. 


A hrst simulation was performed to test all the controllers 
in the task of reference tracking. Zero initial conditions were 
assumed. A reference signal of length 5000 samples (corre¬ 
sponding to 100 s) was used, dehned as a random sequence of 
step signals with amplitudes in the interval [— 7r,7r], hltered by 
a second-order hlter with a cutoff frequency of 10 rad/s. This 
hlter was inserted in order to ensure not too high variations. 
The outputs were corrupted by random uniform noises with 
amplitude 0.02 rad. In Figure]^ the angular positions of the 
closed-loop system with the hrst controller are compared with 
the references for the hrst 20 s of this simulation. Note that 
the two position references were chosen quite similar to each 
other (but not equal) in order to allow the manipulator to reach 
in a simple way any position in its range. A second simulation 
was performed to test the controllers in the task of disturbance 
attenuation. Zero initial conditions and a zero reference were 
assumed. An output disturbance signal of length 1000 samples 
(corresponding to 20 s) was considered, dehned as a sequence 
of two steps (one for each output channel) of amplitude 1 rad, 
hltered by a second-order hlter with a cutoff frequency of 10 
rad/s. The outputs were also corrupted by random uniform 
noises with amplitude 0.02 rad. In Figure the angular 
positions of the closed-loop system with the hrst controller 
are shown, together with the disturbance signals. 


Then, a Monte Carlo (MC) simulation was carried out, where 
this procedure (data generation, control design, reference 
tracking test) was repeated 200 times. For each trial, the 
tracking performance was evaluated by means of the Root 
Mean Square tracking errors, dehned as 



where j is the Ah component of the reference signal and 
yi^t is the Ah component of the controlled system output. 
The average errors RMSi obtained in the MC simulation are 
reported in Table |I^ From these results, it can be concluded 
that the designed control systems are quite effective, showing 
a fast and precise tracking, and a signihcant disturbance 
attenuation capability. In comparison with the two-step method 
of a, the proposed approach is simpler, since a polynomial 
model of the form ([^ has in general a signihcantly simpler 
structure wrt an LPV model (and, in particular, wA a state- 
space LPV model). Moreover, the tracking results obtained 
by the inversion-based controllers are similar (or even slightly 
better) than those obtained by the GS controller, despite the 
fact that this latter uses a stronger information on the system 
^ (i.e., the information that ( [TOl i is a quasi-LPV system). 
The computational times for the control design phase (re¬ 
ferred to a laptop with an i7 3Ghz processor and 16 MB 
RAM) resulted quite low, considering that the set used for 
design consists of 5000 data; 92 s (nonlinear model), 83 s 
(affine model). The control algorithm on-line evaluation times 
resulted also quite low; 2.1e-3 s (nonlinear model), LOe-3 
s (affine model). This shows that these algorithms can be 
effectively implemented on real time processors. 



controller 1 

controller 2 

GS 

RMSi 

0.159 

0.160 

0.167 

RMS 2 

0.114 

0.115 

0.152 


Table II 

Robot Manipulator. Average RMS tracking errors. 



Time [s] 


Figure 4. Robot Manipulator. Continuous (black) line: reference. Dashed 
(red) line: closed-loop system output. 


C. Type 1 diabetes 

A model representing a type 1 diabetic patient has been 
considered in this example. The inputs of this model are the 
carbohydrate-based meal input and the insulin input function, 
the output is the blood glucose concentration (glycemic re¬ 
sponse). The model state equations are the following; 
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Figure 5. Robot Manipulator. Continuous (blue) line: disturbance. Dashed 
(red) line: closed-loop system output. 


= -[pi+ii{t)]yi.t)+piGb +^w{t) 

= -P2V{t) +P3[I{t) - h] 

= - kj{t) (12) 

= + ^u{t) 

= k2lll{t) - {kd + ka)l2{t) 

where y{t) is the blood glucose concentration (the system 
output), I{t) is the blood insulin concentration, r]{t) is the 
insulin concentration in a remote compartment, Vg is the 
volume distribution, w{t) is the carbohydrate-based meal input 
(the system unmeasured input), Ii{t) is the subcutaneous 
insulin mass in the injection depot, l 2 {t) is the subcutaneous 
insulin mass proximal to plasma and u{t) is the injected insulin 
rate (the system measured input); pi, p 2 , ps are individual 
subject parameters, Vd is the plasma distribution volume, ^ 21 , 
ka, kd, and ke are insulin pharmacokinetic parameters, is 
the basal blood insulin concentration and Gh is the basal blood 
glucose concentration. 

The first two equations of ( [T^ , describing the glucose dy¬ 
namics, have been taken from the Bergman model, IT]; the 
last three equations of ( [T2| , describing the insulin kinetics, 
have been taken from the Shimoda model, a. The following 
parameter values have been assumed: pi = 0.031 min~^, 
P 2 = 0.012 p^ = 9.56e — 6 mL/pU, Vg = 

1.45 dL/kg, Vd = 0.2 mL/kg, Vj = 5e — 3 mL, k 2 i = 
0.0166 min~^, ka = 0.0133 min~^, kd — 0.0033 min~^, 
ke = 0.3 h — 0 pU/mL and Gb — 180 mg/dL. In 

this simulated example, the model ( [T^ represents the unknown 
“true” patient metabolic system to control. 

It must be remarked that the model ([T2| is not the most 
recent that can be found in the literature and may also be 
not sufficiently adequate to describe a real diabetes patient. 
However, the aim of this numerical example is to test the 
proposed control algorithm on a non trivial nonlinear system 
and thus the particular choice of the model used as the “true” 
system is not relevant. 


A simulation of the patient system ( [T^ was performed, where 
the insulin input was taken from a set of experimental data, 
measured on a real patient. The meal input was simulated as 
a superposition (with positive coefficients) of exponentially 
decaying signals Wj{t) j = 1, 2,..., where each contribution 
Wj{t) represents a single meal. These signals are of the form 


Wg{t) = 


0, t < tj 

t>tj 


( 13 ) 


where tj is the time at which the patient started to eat. 
The times tj were realistically chosen in order to have an 
insulin injection a few minutes before a meal. A negative 
term of the form ( [T3| ) were also added to the meal input in 
order to reproduce the effects of an external input yielding a 
decrease of the output (e.g. a physical activity). The output 
signal (the blood glucose concentration) resulting from this 
simulation was corrupted by a white noise, having a noise-to- 
signal standard deviation ratio of 3%. 

A set of L = 4800 data (corresponding to 10 days) was 
collected from this simulation, using a sampling time Tg = 3 
min'. 

'D = {Wt, yt}t=_47gg 


where Ut = u{Tgt) are the measurements of the insulin input 
and yt = yiTgt) are the measurements of the output. Note 
that, as it happens in most real situations, the meal input was 
not measured. 


A nonlinear controller was designed following the approach 
described in Sections and |I^ This controller was applied 
to the diabetes system ( fTTl i. 

Three simulations of the patient system ( |T^ with duration 10 
days were performed, using a meal input signal different from 
that used to generate the design data T). The insulin signal 
was generated as follows: 

• first simulation: zero insulin; 

• second simulation: insulin injected by the patient on the 
basis of his/hers experience; 

• third simulation: insulin signal computed by the designed 
controller. 

In the simulations, the output signal was corrupted by a white 
noise, with a noise-to-signal standard deviation ratio of 3%. 
The obtained results can be commented as follows: With no 
insulin, the glucose concentration becomes very large, leading 
to serious health problems of the patient. When the amount 
of injected insulin is decided by the patient, the glucose 
concentration is somewhat regulated but it may reach large 
values, which may worsen the patient health conditions (see 
Figure |^. When the amount of injected insulin is decided by 
the controller, the glucose concentration is always kept within 
the interval [80,180] mg/dL which, in diabetes treatment 
medicine, is commonly considered a safe interval (see Figure 
1 ^. 
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Figure 6. Blue lines: patient glucose concentration; red lines: safety bounds. 
Above: result with insulin signal decided by the patient; below: result with 
insulin signal generated by the controller. 
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